"""
FC Exclusion Decomposition: Why does broad financial center exclusion
kill the KAOPEN interactions in the gravity bilateral followup?

Hypothesis (a): Result driven by a handful of financial hubs.
Hypothesis (b): Removing high-KAOPEN countries compresses variation,
                mechanically destroying interaction power.

Tests:
1. KAOPEN_j distribution across exclusion sets
2. Intermediate exclusions (SGP+HKG, CHE+NLD+BEL, GBR)
3. Broad exclusion sample KAOPEN range/SD
4. FC countries as reporters only vs partners only
5. Artificial KAOPEN compression test (trim top KAOPEN in full sample)
"""

import pandas as pd
import numpy as np
from pathlib import Path
import sys
from io import StringIO

sys.path.insert(0, str(Path("/mnt/c/demographics_capital_flows/gravity_bilateral_followup")))
from src.model import PanelGLS

BASE_DIR = Path("/mnt/c/demographics_capital_flows/gravity_bilateral_followup")
PROCESSED_DIR = BASE_DIR / "data" / "processed"
OUTPUT_DIR = BASE_DIR / "output" / "tables"

FINANCIAL_CENTERS_NARROW = {
    'LUX', 'IRL', 'CYM', 'BMU', 'BHS', 'PAN', 'VGB', 'BHR', 'MUS', 'MLT', 'CYP'
}
FINANCIAL_CENTERS_BROAD = FINANCIAL_CENTERS_NARROW | {
    'HKG', 'SGP', 'CHE', 'NLD', 'BEL', 'GBR'
}
# The 6 countries added in broad vs narrow
BROAD_ADDITIONS = {'HKG', 'SGP', 'CHE', 'NLD', 'BEL', 'GBR'}


def estimate_model_2c(df, label):
    """Estimate Model 2c and return key results."""
    gravity_vars = ['log_dist', 'contiguity', 'common_lang_official',
                    'colonial_ties', 'log_gdp_product']
    demo_vars = ['dZ_1', 'dZ_2', 'dZ_3']
    kaopen_vars = ['kaopen_j', 'dZ_1_x_kaopen_j', 'dZ_2_x_kaopen_j', 'dZ_3_x_kaopen_j']
    dep_var = 'log_portfolio_total'

    years = sorted(df['year'].dropna().unique())
    yr_cols = [f'yr_{int(y)}' for y in years[1:]]
    yr_cols = [c for c in yr_cols if c in df.columns]

    regressors = gravity_vars + demo_vars + kaopen_vars
    all_vars = regressors + yr_cols

    est = df.dropna(subset=[dep_var] + all_vars + ['pair_id', 'year']).copy()
    if len(est) < 100:
        return None

    y = est[dep_var].values
    X = est[all_vars].values

    gls = PanelGLS()
    gls.fit(y, X, est['pair_id'].values, est['year'].values.astype(int))

    # Extract key coefficients
    result = {
        'label': label,
        'n_obs': gls.n_obs,
        'n_pairs': gls.n_countries,
        'r_squared': gls.r_squared,
    }
    for i, v in enumerate(regressors):
        result[f'{v}_coef'] = gls.beta[i]
        result[f'{v}_se'] = gls.se[i]
        result[f'{v}_pval'] = gls.pvalues[i]

    # Also return KAOPEN stats from the estimation sample
    result['kaopen_j_mean'] = est['kaopen_j'].mean()
    result['kaopen_j_sd'] = est['kaopen_j'].std()
    result['kaopen_j_min'] = est['kaopen_j'].min()
    result['kaopen_j_max'] = est['kaopen_j'].max()
    result['kaopen_j_p10'] = est['kaopen_j'].quantile(0.10)
    result['kaopen_j_p90'] = est['kaopen_j'].quantile(0.90)
    result['kaopen_j_iqr'] = est['kaopen_j'].quantile(0.75) - est['kaopen_j'].quantile(0.25)
    result['n_unique_partners'] = est['partner'].nunique()
    result['n_high_kaopen_partners'] = est.loc[est['kaopen_j'] > 2.0, 'partner'].nunique()
    result['pct_obs_high_kaopen'] = (est['kaopen_j'] > 2.0).mean() * 100

    # Interaction term variation
    result['interaction_sd'] = est['dZ_1_x_kaopen_j'].std()

    return result


def main():
    print("=" * 70)
    print("FC EXCLUSION DECOMPOSITION")
    print("=" * 70)

    df = pd.read_csv(PROCESSED_DIR / "bilateral_panel.csv")
    print(f"Full panel: {len(df):,} obs")

    # Create year dummies
    years = sorted(df['year'].dropna().unique())
    for y in years[1:]:
        df[f'yr_{int(y)}'] = (df['year'] == y).astype(int)

    all_results = []

    # ===================================================================
    # PART 1: Replicate baselines + intermediate exclusions
    # ===================================================================
    print("\n" + "=" * 70)
    print("PART 1: MODEL ESTIMATES ACROSS EXCLUSION SETS")
    print("=" * 70)

    exclusion_sets = {
        '(1) Full sample': set(),
        '(2) Narrow FC excl': FINANCIAL_CENTERS_NARROW,
        '(3) +SGP+HKG only': FINANCIAL_CENTERS_NARROW | {'SGP', 'HKG'},
        '(4) +CHE+NLD+BEL only': FINANCIAL_CENTERS_NARROW | {'CHE', 'NLD', 'BEL'},
        '(5) +GBR only': FINANCIAL_CENTERS_NARROW | {'GBR'},
        '(6) +SGP+HKG+GBR': FINANCIAL_CENTERS_NARROW | {'SGP', 'HKG', 'GBR'},
        '(7) +CHE+NLD+BEL+GBR': FINANCIAL_CENTERS_NARROW | {'CHE', 'NLD', 'BEL', 'GBR'},
        '(8) Broad FC excl (all 6)': FINANCIAL_CENTERS_BROAD,
    }

    for label, fc_set in exclusion_sets.items():
        if fc_set:
            mask = ~(df['reporter'].isin(fc_set) | df['partner'].isin(fc_set))
            df_sub = df[mask].copy()
        else:
            df_sub = df.copy()

        n_dropped = len(df) - len(df_sub)
        print(f"\n--- {label} ---")
        print(f"  Excluded: {sorted(fc_set) if fc_set else 'none'}")
        print(f"  Dropped: {n_dropped:,} obs ({n_dropped/len(df)*100:.1f}%)")

        res = estimate_model_2c(df_sub, label)
        if res is not None:
            all_results.append(res)
            sig = '***' if res['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if res['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if res['dZ_1_x_kaopen_j_pval'] < 0.1 else 'NS'
            print(f"  dZ_1 x KAOPEN_j = {res['dZ_1_x_kaopen_j_coef']:.3f} (p={res['dZ_1_x_kaopen_j_pval']:.4f}) {sig}")
            print(f"  KAOPEN_j SD={res['kaopen_j_sd']:.3f}, IQR={res['kaopen_j_iqr']:.3f}, "
                  f"interaction SD={res['interaction_sd']:.4f}")
            print(f"  High-KAOPEN partners (>2.0): {res['n_high_kaopen_partners']}, "
                  f"{res['pct_obs_high_kaopen']:.1f}% of obs")

    # ===================================================================
    # PART 2: Reporter vs Partner decomposition
    # ===================================================================
    print("\n" + "=" * 70)
    print("PART 2: REPORTER vs PARTNER DECOMPOSITION")
    print("=" * 70)

    # FC countries as REPORTERS only (keep as partners)
    mask_reporter_only = ~df['reporter'].isin(BROAD_ADDITIONS)
    df_excl_reporter = df[mask_reporter_only].copy()
    n_drop = len(df) - len(df_excl_reporter)
    print(f"\n--- Excl broad-6 as REPORTERS only (keep as partners) ---")
    print(f"  Dropped: {n_drop:,} obs ({n_drop/len(df)*100:.1f}%)")
    res = estimate_model_2c(df_excl_reporter, 'Excl broad-6 as reporters')
    if res is not None:
        all_results.append(res)
        sig = '***' if res['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if res['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if res['dZ_1_x_kaopen_j_pval'] < 0.1 else 'NS'
        print(f"  dZ_1 x KAOPEN_j = {res['dZ_1_x_kaopen_j_coef']:.3f} (p={res['dZ_1_x_kaopen_j_pval']:.4f}) {sig}")

    # FC countries as PARTNERS only (keep as reporters)
    mask_partner_only = ~df['partner'].isin(BROAD_ADDITIONS)
    df_excl_partner = df[mask_partner_only].copy()
    n_drop = len(df) - len(df_excl_partner)
    print(f"\n--- Excl broad-6 as PARTNERS only (keep as reporters) ---")
    print(f"  Dropped: {n_drop:,} obs ({n_drop/len(df)*100:.1f}%)")
    res = estimate_model_2c(df_excl_partner, 'Excl broad-6 as partners')
    if res is not None:
        all_results.append(res)
        sig = '***' if res['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if res['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if res['dZ_1_x_kaopen_j_pval'] < 0.1 else 'NS'
        print(f"  dZ_1 x KAOPEN_j = {res['dZ_1_x_kaopen_j_coef']:.3f} (p={res['dZ_1_x_kaopen_j_pval']:.4f}) {sig}")

    # ONLY FC countries as reporters (rest of world as partners)
    mask_fc_reporter = df['reporter'].isin(BROAD_ADDITIONS) & ~df['partner'].isin(FINANCIAL_CENTERS_NARROW)
    df_fc_reporter = df[mask_fc_reporter].copy()
    print(f"\n--- ONLY broad-6 as REPORTERS (non-narrow-FC partners) ---")
    print(f"  Sample: {len(df_fc_reporter):,} obs")
    res = estimate_model_2c(df_fc_reporter, 'Only broad-6 as reporters')
    if res is not None:
        all_results.append(res)
        sig = '***' if res['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if res['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if res['dZ_1_x_kaopen_j_pval'] < 0.1 else 'NS'
        print(f"  dZ_1 x KAOPEN_j = {res['dZ_1_x_kaopen_j_coef']:.3f} (p={res['dZ_1_x_kaopen_j_pval']:.4f}) {sig}")
    else:
        print("  Insufficient observations")

    # ONLY FC countries as partners (rest of world as reporters)
    mask_fc_partner = df['partner'].isin(BROAD_ADDITIONS) & ~df['reporter'].isin(FINANCIAL_CENTERS_NARROW)
    df_fc_partner = df[mask_fc_partner].copy()
    print(f"\n--- ONLY broad-6 as PARTNERS (non-narrow-FC reporters) ---")
    print(f"  Sample: {len(df_fc_partner):,} obs")
    res = estimate_model_2c(df_fc_partner, 'Only broad-6 as partners')
    if res is not None:
        all_results.append(res)
        sig = '***' if res['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if res['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if res['dZ_1_x_kaopen_j_pval'] < 0.1 else 'NS'
        print(f"  dZ_1 x KAOPEN_j = {res['dZ_1_x_kaopen_j_coef']:.3f} (p={res['dZ_1_x_kaopen_j_pval']:.4f}) {sig}")
    else:
        print("  Insufficient observations")

    # ===================================================================
    # PART 3: KAOPEN variation analysis
    # ===================================================================
    print("\n" + "=" * 70)
    print("PART 3: KAOPEN_j DISTRIBUTION ANALYSIS")
    print("=" * 70)

    dep_var = 'log_portfolio_total'
    est_vars = ['log_dist', 'contiguity', 'common_lang_official', 'colonial_ties',
                'log_gdp_product', 'dZ_1', 'dZ_2', 'dZ_3', 'kaopen_j',
                'dZ_1_x_kaopen_j', 'dZ_2_x_kaopen_j', 'dZ_3_x_kaopen_j', 'pair_id', 'year']
    est_full = df.dropna(subset=[dep_var] + est_vars).copy()

    for label, fc_set in [('Full sample', set()),
                           ('Narrow FC excl', FINANCIAL_CENTERS_NARROW),
                           ('Broad FC excl', FINANCIAL_CENTERS_BROAD)]:
        if fc_set:
            mask = ~(est_full['reporter'].isin(fc_set) | est_full['partner'].isin(fc_set))
            sub = est_full[mask]
        else:
            sub = est_full

        print(f"\n  {label}:")
        print(f"    N obs: {len(sub):,}")
        print(f"    N unique partners: {sub['partner'].nunique()}")
        k = sub['kaopen_j']
        print(f"    KAOPEN_j: mean={k.mean():.3f}, sd={k.std():.3f}, "
              f"min={k.min():.3f}, max={k.max():.3f}")
        print(f"    KAOPEN_j quantiles: p10={k.quantile(.1):.3f}, p25={k.quantile(.25):.3f}, "
              f"p50={k.quantile(.5):.3f}, p75={k.quantile(.75):.3f}, p90={k.quantile(.9):.3f}")

        # Which partners have KAOPEN > 2.0?
        high_k = sub[sub['kaopen_j'] > 2.0]
        partners_high = high_k['partner'].unique()
        print(f"    Partners with KAOPEN_j > 2.0: {len(partners_high)} "
              f"({', '.join(sorted(partners_high)[:10])}{'...' if len(partners_high)>10 else ''})")
        print(f"    Obs with KAOPEN_j > 2.0: {len(high_k):,} ({len(high_k)/len(sub)*100:.1f}%)")

        # dZ_1 x kaopen_j variation
        ix = sub['dZ_1_x_kaopen_j']
        print(f"    dZ_1 x KAOPEN_j: sd={ix.std():.4f}, range=[{ix.min():.3f}, {ix.max():.3f}]")

    # ===================================================================
    # PART 4: Country-by-country KAOPEN values for the 6 broad additions
    # ===================================================================
    print("\n" + "=" * 70)
    print("PART 4: KAOPEN VALUES OF THE 6 BROAD-ADDITION COUNTRIES")
    print("=" * 70)

    for iso in sorted(BROAD_ADDITIONS):
        sub = est_full[est_full['partner'] == iso]
        if len(sub) > 0:
            k_vals = sub.groupby('year')['kaopen_j'].first()
            print(f"\n  {iso}: {len(sub):,} obs as partner, KAOPEN range [{k_vals.min():.3f}, {k_vals.max():.3f}], "
                  f"mean={k_vals.mean():.3f}")
        else:
            print(f"\n  {iso}: no observations as partner")

    # ===================================================================
    # PART 5: Artificial KAOPEN compression test
    # ===================================================================
    print("\n" + "=" * 70)
    print("PART 5: ARTIFICIAL KAOPEN COMPRESSION TEST")
    print("=" * 70)
    print("  If we keep the full sample but artificially cap KAOPEN_j at the")
    print("  broad-exclusion maximum, does the interaction die?")

    # Get broad-exclusion max KAOPEN
    broad_mask = ~(est_full['reporter'].isin(FINANCIAL_CENTERS_BROAD) |
                   est_full['partner'].isin(FINANCIAL_CENTERS_BROAD))
    broad_sub = est_full[broad_mask]
    broad_kaopen_max = broad_sub['kaopen_j'].max()
    broad_kaopen_p90 = broad_sub['kaopen_j'].quantile(0.90)
    print(f"  Broad exclusion KAOPEN_j max: {broad_kaopen_max:.3f}")

    # Test: cap KAOPEN_j at broad-exclusion max in full sample
    df_capped = df.copy()
    df_capped['kaopen_j'] = df_capped['kaopen_j'].clip(upper=broad_kaopen_max)
    df_capped['dZ_1_x_kaopen_j'] = df_capped['dZ_1'] * df_capped['kaopen_j']
    df_capped['dZ_2_x_kaopen_j'] = df_capped['dZ_2'] * df_capped['kaopen_j']
    df_capped['dZ_3_x_kaopen_j'] = df_capped['dZ_3'] * df_capped['kaopen_j']
    res = estimate_model_2c(df_capped, 'Full sample, KAOPEN capped at broad max')
    if res is not None:
        all_results.append(res)
        sig = '***' if res['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if res['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if res['dZ_1_x_kaopen_j_pval'] < 0.1 else 'NS'
        print(f"  dZ_1 x KAOPEN_j (capped) = {res['dZ_1_x_kaopen_j_coef']:.3f} (p={res['dZ_1_x_kaopen_j_pval']:.4f}) {sig}")

    # ===================================================================
    # PART 6: Leave-one-out from broad additions
    # ===================================================================
    print("\n" + "=" * 70)
    print("PART 6: LEAVE-ONE-OUT (drop each of 6 broad countries individually)")
    print("=" * 70)

    for iso in sorted(BROAD_ADDITIONS):
        fc_set = FINANCIAL_CENTERS_NARROW | {iso}
        mask = ~(df['reporter'].isin(fc_set) | df['partner'].isin(fc_set))
        df_sub = df[mask].copy()
        res = estimate_model_2c(df_sub, f'Narrow + excl {iso}')
        if res is not None:
            all_results.append(res)
            sig = '***' if res['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if res['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if res['dZ_1_x_kaopen_j_pval'] < 0.1 else 'NS'
            print(f"  Excl {iso}: dZ_1 x KAOPEN_j = {res['dZ_1_x_kaopen_j_coef']:.3f} "
                  f"(p={res['dZ_1_x_kaopen_j_pval']:.4f}) {sig}  "
                  f"[N={res['n_obs']:,}, KAOPEN SD={res['kaopen_j_sd']:.3f}]")

    # ===================================================================
    # BUILD OUTPUT TABLE
    # ===================================================================
    results_df = pd.DataFrame(all_results)

    # Save raw results
    results_df.to_csv(OUTPUT_DIR / "fc_exclusion_decomposition_raw.csv", index=False)

    # Build markdown report
    md = StringIO()
    md.write("# Financial Center Exclusion Decomposition\n\n")
    md.write("**Question**: Why does broad FC exclusion kill the KAOPEN interaction?\n\n")
    md.write("- Full sample: dZ_1 x KAOPEN_j = 0.967 (p=0.014)\n")
    md.write("- Narrow FC exclusion: dZ_1 x KAOPEN_j = 0.826 (p=0.049)\n")
    md.write("- Broad FC exclusion: dZ_1 x KAOPEN_j = 0.386 (p=0.414)\n\n")

    # Table 1: Intermediate exclusions
    md.write("## Table 1: Intermediate Exclusions (Model 2c)\n\n")
    md.write("| Exclusion | N obs | dZ1xKAOPEN coef | p-value | KAOPEN_j SD | Interaction SD | High-KAOPEN partners |\n")
    md.write("|-----------|------:|----------------:|--------:|------------:|---------------:|---------------------:|\n")
    for r in all_results:
        label = r['label']
        if label.startswith('(') or label == 'Full sample, KAOPEN capped at broad max':
            sig = '***' if r['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if r['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if r['dZ_1_x_kaopen_j_pval'] < 0.1 else ''
            md.write(f"| {label} | {r['n_obs']:,} | {r['dZ_1_x_kaopen_j_coef']:.3f}{sig} | "
                     f"{r['dZ_1_x_kaopen_j_pval']:.4f} | {r['kaopen_j_sd']:.3f} | "
                     f"{r['interaction_sd']:.4f} | {r['n_high_kaopen_partners']} |\n")
    md.write("\n")

    # Table 2: Reporter vs Partner
    md.write("## Table 2: Reporter vs Partner Decomposition\n\n")
    md.write("| Sample | N obs | dZ1xKAOPEN coef | p-value | KAOPEN_j SD |\n")
    md.write("|--------|------:|----------------:|--------:|------------:|\n")
    for r in all_results:
        label = r['label']
        if 'reporter' in label.lower() or 'partner' in label.lower():
            sig = '***' if r['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if r['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if r['dZ_1_x_kaopen_j_pval'] < 0.1 else ''
            md.write(f"| {label} | {r['n_obs']:,} | {r['dZ_1_x_kaopen_j_coef']:.3f}{sig} | "
                     f"{r['dZ_1_x_kaopen_j_pval']:.4f} | {r['kaopen_j_sd']:.3f} |\n")
    md.write("\n")

    # Table 3: Leave-one-out
    md.write("## Table 3: Leave-One-Out from Broad Additions\n\n")
    md.write("| Excluded country | N obs | dZ1xKAOPEN coef | p-value | KAOPEN_j SD |\n")
    md.write("|-----------------|------:|----------------:|--------:|------------:|\n")
    for r in all_results:
        label = r['label']
        if label.startswith('Narrow + excl'):
            iso = label.split('excl ')[-1]
            sig = '***' if r['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if r['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if r['dZ_1_x_kaopen_j_pval'] < 0.1 else ''
            md.write(f"| {iso} | {r['n_obs']:,} | {r['dZ_1_x_kaopen_j_coef']:.3f}{sig} | "
                     f"{r['dZ_1_x_kaopen_j_pval']:.4f} | {r['kaopen_j_sd']:.3f} |\n")
    md.write("\n")

    # Table 4: KAOPEN compression test
    md.write("## Table 4: Artificial KAOPEN Compression Test\n\n")
    for r in all_results:
        if r['label'] == 'Full sample, KAOPEN capped at broad max':
            sig = '***' if r['dZ_1_x_kaopen_j_pval'] < 0.01 else '**' if r['dZ_1_x_kaopen_j_pval'] < 0.05 else '*' if r['dZ_1_x_kaopen_j_pval'] < 0.1 else ''
            md.write(f"Full sample with KAOPEN_j capped at {broad_kaopen_max:.3f} (broad-exclusion max):\n")
            md.write(f"- dZ_1 x KAOPEN_j = {r['dZ_1_x_kaopen_j_coef']:.3f} (p={r['dZ_1_x_kaopen_j_pval']:.4f}) {sig}\n")
            md.write(f"- N obs: {r['n_obs']:,}\n\n")

    # Diagnosis
    md.write("## Diagnosis\n\n")

    # Collect key numbers for diagnosis
    full_res = next((r for r in all_results if r['label'] == '(1) Full sample'), None)
    narrow_res = next((r for r in all_results if r['label'] == '(2) Narrow FC excl'), None)
    broad_res = next((r for r in all_results if r['label'] == '(8) Broad FC excl (all 6)'), None)
    capped_res = next((r for r in all_results if r['label'] == 'Full sample, KAOPEN capped at broad max'), None)

    if full_res and broad_res:
        sd_drop = (full_res['kaopen_j_sd'] - broad_res['kaopen_j_sd']) / full_res['kaopen_j_sd'] * 100
        ix_sd_drop = (full_res['interaction_sd'] - broad_res['interaction_sd']) / full_res['interaction_sd'] * 100
        partner_drop = full_res['n_high_kaopen_partners'] - broad_res['n_high_kaopen_partners']

        md.write(f"**KAOPEN_j variation lost**: SD drops {sd_drop:.1f}% (from {full_res['kaopen_j_sd']:.3f} to {broad_res['kaopen_j_sd']:.3f})\n\n")
        md.write(f"**Interaction term variation lost**: SD drops {ix_sd_drop:.1f}% (from {full_res['interaction_sd']:.4f} to {broad_res['interaction_sd']:.4f})\n\n")
        md.write(f"**High-KAOPEN partners lost**: {partner_drop} countries with KAOPEN_j > 2.0\n\n")

    if capped_res:
        md.write(f"**Compression test**: When KAOPEN_j is capped at {broad_kaopen_max:.3f} in the full sample, ")
        if capped_res['dZ_1_x_kaopen_j_pval'] < 0.05:
            md.write(f"the interaction SURVIVES (p={capped_res['dZ_1_x_kaopen_j_pval']:.4f}). ")
            md.write("This means removing the countries themselves (not just their KAOPEN variation) is what kills the result.\n\n")
            md.write("**Conclusion: Hypothesis (a) -- the result IS driven by flows to/from financial hubs.**\n")
        elif capped_res['dZ_1_x_kaopen_j_pval'] > 0.10:
            md.write(f"the interaction ALSO DIES (p={capped_res['dZ_1_x_kaopen_j_pval']:.4f}). ")
            md.write("This means the interaction needs the high-KAOPEN variation that these countries provide.\n\n")
            md.write("**Conclusion: Hypothesis (b) -- removing high-KAOPEN observations mechanically destroys identification.**\n")
        else:
            md.write(f"the interaction is marginal (p={capped_res['dZ_1_x_kaopen_j_pval']:.4f}). ")
            md.write("The truth is likely a mix of both hypotheses.\n\n")

    # Check reporter vs partner
    rpt_res = next((r for r in all_results if r['label'] == 'Excl broad-6 as reporters'), None)
    prt_res = next((r for r in all_results if r['label'] == 'Excl broad-6 as partners'), None)
    if rpt_res and prt_res:
        md.write(f"\n**Reporter vs Partner**: ")
        if rpt_res['dZ_1_x_kaopen_j_pval'] < 0.05 and prt_res['dZ_1_x_kaopen_j_pval'] > 0.10:
            md.write("Effect survives when excluding as reporters but dies when excluding as partners. ")
            md.write("The interaction is driven by these countries as DESTINATIONS (high-KAOPEN partners attract more demographic-driven flows).\n")
        elif prt_res['dZ_1_x_kaopen_j_pval'] < 0.05 and rpt_res['dZ_1_x_kaopen_j_pval'] > 0.10:
            md.write("Effect survives when excluding as partners but dies when excluding as reporters. ")
            md.write("The interaction is driven by these countries as SOURCES (financial hubs send more demographic-driven flows).\n")
        elif rpt_res['dZ_1_x_kaopen_j_pval'] < 0.05 and prt_res['dZ_1_x_kaopen_j_pval'] < 0.05:
            md.write("Effect survives both exclusions. The interaction is robust; the full broad exclusion simply removes too many obs.\n")
        else:
            md.write(f"Reporter-excl p={rpt_res['dZ_1_x_kaopen_j_pval']:.4f}, "
                     f"Partner-excl p={prt_res['dZ_1_x_kaopen_j_pval']:.4f}. "
                     f"Both weaken, consistent with these countries being important on both sides.\n")

    outfile = OUTPUT_DIR / "fc_exclusion_decomposition.md"
    outfile.write_text(md.getvalue())
    print(f"\nSaved: {outfile}")
    print(f"Saved: {OUTPUT_DIR / 'fc_exclusion_decomposition_raw.csv'}")


if __name__ == "__main__":
    main()
